Three-dimensional dendrite tip morphology at low undercooling 
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We investigate the three-dimensional morphology of the dendrite tip using the phase-field method. 
We find that, for low undercoolings, this morphology is ostensibly independent of anisotropy strength 
except for a localized shape distortion near the tip that only affects the value of the tip radius p 
(which is crudely approximated by p ~ (1 — a)pi v where pi v is the Ivantsov tip radius of an 
isothermal paraboloid with the same tip velocity and a is the stiffness anisotropy). The universal 
tip shape, which excludes this distortion, is well fitted by the form z = —r 2 /2 + A^r 4, cos 40 where 
\z\ is the distance from the tip and all lengths are scaled by pi v . This fit yields A4 in the range 
0.004 — 0.005 in good quantitative agreement with the existing tip morphology measurements in 
succinonitrile [LaCombe et al, Phys. Rev. E 52, 2778 (1995)], which are reanalyzed here and found 
to be consistent with a single cos 40 mode non-axisymmetric deviation from a paraboloid. Moreover, 
the fin shape away from the tip is well fitted by the power law z = — a\x\ 5 ^ 3 with a ta 0.68. Finally, 
the characterization of the operating state of the dendrite tip is revisited in the light of these results. 
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I. INTRODUCTION 

The shape of crystal dendrites was first suggested by 
Papapetrou (!]] to be parabolic, and slightly more than a 
decade later Ivantsov demonstrated [Q] that a parabola 
(paraboloid) is an exact solution of the steady-state 
growth equations in two (three) dimensions when cap- 
illary effects are entirely neglected and the interface is 
isothermal. He derived the well-known relationship 



A = P Iv exp(Pr w ) 



ds 



exp( 



(1) 



Plr 



between the Peclet number Pj v = pj v V/(2D) and the 
dimensionless undercooling A = (Tm ~ T^)/ (L/c p ) for a 
paraboloid of tip velocity V and tip radius pr v where D 
is the thermal diffusivity, Tm is the melting temperature, 
Too is the initial temperature of the undercooled liquid, 
L is the latent heat of melting, and c p is the specific heat 
at constant pressure. 

In more recent history, the development of solvability 
theory has led to the additional and crucial under- 
standing that the anisotropic surface energy acts as a sin- 
gular perturbation that uniquely selects the tip velocity, 
and also alters the entire dendrite shape . According 
to this theory, the scaling parameter a* = 2Ddo/ pj v V 
(where do, defined below, is the capillary length) ap- 
proaches a constant at low undercooling that only de- 
pends on the anisotropy strength (denoted here by 64 for 
a crystal with an underlying cubic symmetry) , and scales 

as a* ~ e 7 J 4 in the limit of vanishingly small anisotropy. 
In this same limit, Ben Amar and Brener || have pre- 
dicted that capillary effects lead to a universal four-fold 
deviation from a paraboloid of the form 



where A4 — 1/88 is independent of anisotropy strength 
and (r, <f>) are the polar coordinates in the plane normal 
to the growth axis z with all lengths scaled by pj v . The 
improved prediction A4 = 1/96 has been obtained in a 
subsequent analysis 

To characterize the shape further behind the tip, 
Brener remarked that the cross-sectional shape of nee- 
dle crystal dendrites perpendicular to the growth axis 
can be assumed to evolve with increasing distance \z\ 
from the tip as a two-dimensional growth shape in time 
(i = |z|/V). The area of this shape increases linearly in 
time, and thus \z\, as 



S(z) = 2np Iv \z\ 



(3) 



where pj v is defined by Eq. [l] with V being the ac- 
tual tip velocity. This area corresponds to the well- 
known self-similarity solution of a growing circle in two 
dimensions and is also the exact cross-sectional area of 
the three-dimensional Ivantsov paraboloid of revolution. 
The above mapping is presumed to become justified far 
enough away from the tip where the component of the 
heat flux along z can be neglected. Using this mapping 
and the results of a previous analysis of two-dimensional 
anisotropic Laplacian growth at constant flux |^| , he pre- 
dicted that the width of the four dendrite fins should 
increase as the 3/5 power of the distance behind the tip 
, yielding the expression 



1/3 



1 5/3 



(4) 



z = — — + Air 4 cos4</>, 



(2) 



where x is the interface coordinate normal to the z-axis 
in a (010) section of the tip and crj ( <t *) i s the tip scal- 
ing parameter in two (three) dimensions. Moreover, the 
ratio el / a * becomes independent of anisotropy for weak 
anisotropy with (<7| /er*) 1 / 3 ~ 1. 
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Our goal in this paper is to use the phase-field method 
to obtain a characterization of the dendrite tip morphol- 
ogy that is sufficiently detailed and accurate to test the 
above theoretical predictions and to make a critical com- 
parison with existing benchmark experiments PJTc|] . The 
present simulations are based on a novel adaptive-step 
diffusion Monte Carlo method [ pd|Jl2| that provides an 
efficient treatment of the large scale diffusion field away 
from the growing structure. Thus this method allows us 
to investigate a relatively low undercooling range A ~ 0.1 
(Pj v ~ 0.03), in contrast to previous simulations that 
were limited to Peclet numbers of order unity |13|. As 
it turns out, and this is one of the main findings in this 
paper, the tip morphology in this Peclet number range 
is already ostensibly independent of anisotropy strength 
and undercooling, except for a localized shape distortion 
at the tip that is not experimentally relevant. We are 
therefore able to compare meaningfully this morphology 
to existing detailed dendrite shape measurements in suc- 
cinonitrile P] and xenon jl0| , even though these measure- 
ments were carried out at even smaller A. 

In the next section, we review the basic equations of 
the symmetric model of dendritic growth in three dimen- 
sions and briefly summarize our numerical methods. The 
numerical results are then presented in section III and 
discussed in section IV. Finally, a summary and conclu- 
sions are presented in section V. 



II. BASIC EQUATIONS 

We study the standard symmetric model of solidifi- 
cation in a pure undercooled melt that assumes equal 
thermal diffusivities in the solid and liquid phases. The 
basic equations of this model are given by 



d t u = DV 2 u, 
v n = D(d n u\ 



<9 n w| ; ) . 



-do y~] [a(n) + dg.a(h) 



(5) 
(6) 

(7) 



i=i 



where following common notation, u = (T — T m )/(L/c p ) 
is the scaled temperature field that is zero in equilibrium 
and equal to — A in the liquid far from the interface, d n u\i 
(d n u\ s ) is the normal gradient of u on the liquid (solid) 
side of the interface, v n is the normal velocity of the inter- 
face, 9i are the local angles between the normal h to the 
interface and the two local principal directions on the in- 
terface, Ki are the principal curvatures, do = 7oTmc p /Z/ 2 
is the microscopic capillary length, and 7(n) = 7oa(n) is 
the anisotropic surface energy where 



a(n) = (1 — 364) 



4e 4 



1 -3e 4 



(8) 



for a material with an underlying cubic symmetry, with 
the x, y, and z axes chosen parallel to the [100] direc- 
tions. In writing down the interface condition (0), we 




FIG. 1. Snapshot of a three-dimensional simulated den- 
drite with A = 0.1 and e 4 = 0.025. 



have purposely neglected the effect of interface kinetics 
that is believed to be negligible for the low undercooling 
range where benchmark shape measurements have been 
carried out PjiQl. Moreover, theoretical shape predic- 
tions to date |5|-|7|] have neglected this effect. 

We use two different numerical methods. The first, and 
main one used here, is the phase-field approach JTi 15 
that allows us to study the dendrite shape for the full cu- 
bic form of the surface energy defined by Eq. ^|. In partic- 
ular, we simulate the time-dependent evolution of a single 
dendrite arm until steady-state growth is reached. Since 
we are primarily interested in a low undercooling regime, 
we use the results of the asymptotic analysis of Karma 
and Rappel Ea] in order to choose the interface width 
in the phase-field model about an order of magnitude 
smaller than p, and thus much larger than do for small 
A. This analysis also allows us to choose the model pa- 
rameters so as to make kinetic effects vanishingly small. 
Moreover, we use a recently developed Monte Carlo al- 
gorithm to integrate efficiently the diffusion equation in 
the liquid away from the growing structure. In this ap- 
proach, space is divided into two regions. The first region 
consists of the solid plus a thin liquid layer surrounding 
the interface where the deterministic phase-field equa- 
tions are solved on a cubic lattice with the same choice 
of computational parameters as in Ref. jL3|. The sec- 
ond region is the rest of the liquid where the diffusion 
equation is solved using an ensemble of random walk- 
ers that take progressively larger steps with increasing 
distance away from the interface. The algorithm used 
to interface the deterministic and stochastic solutions of 
the diffusion equation in these two regions and to update 
the walkers has been summarized in Ref. [[llj where it 
was used to study the early stage of dendritic evolution. 
It has also been exposed in more details in Ref. [0| to- 
gether with the results of numerical tests and need not 
be re-described here. 
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We exploit the cubic symmetry to reduce simulation 
time by only integrating a part of a single dendrite arm 
(the domain x>0, 0<y<x and z > x), taking 
advantage of the symmetry planes defined by x = 0, y = 
0, x = y, and x — z. The whole dendrite can then be 
reconstructed by successive reflections at these planes, 
and the result of one of our simulations is shown in Fig. 
[j]. These simulations were performed on regular lattices 
of size 240 x 240 x 800 and each took about 200 hours 
of CPU time on a 525 MHz DEC-8400 computer. Note 
that almost no sidebranches can be discerned, although 
the length of the dendrite arms is more than 40 times the 
tip radius of curvature. For the analysis of the steady- 
state shapes, we used only the part of the dendrite grown 
at a constant tip velocity, which corresponds to about 
one third of the arm length in Fig. [lj. To check whether 
the shape is well converged, we performed one run in 
which the dendrite arm growing along the positive z- 
direction was "cut off" sufficiently far behind the tip, 
and its further evolution was simulated in a smaller box 
moving with the dendrite tip for several diffusion times 
(D/V 2 ). No significant change of the shape was observed. 
In addition, this run (for A = 0.1 and £4 = 0.025) yielded 
data which extend considerably farther behind the tip 
than for the other parameters, well to within the range 
where sidebranches are observed in experiments. 

The second numerical method is the standard bound- 
ary integral method (r| that can be used to solve di- 
rectly the sharp-interface steady-state growth equations 
with the simplified axisymmctric form of surface energy 



a(h) = (1 — 3e4) 



1 - 



4e 4 



l-3e 4 



(cos 4 1 



• sin 4 9) 



> (9) 



which is obtained by averaging the full cubic form (^) 
over the polar angle <f> in the x — y plane; 9 is the angle 
between the local normal to the solid-liquid interface and 
the z-axis. This method only describes axisymmetric tip 
shapes and is thus only used here as an additional ba- 
sis of comparison with phase-field results regarding the 
anisotropy-dependent shape distortion in the tip neigh- 
borhood and the selected tip operating state. 



III. NUMERICAL RESULTS 

A. Equal cross-sectional area shape 

The cross-sectional area of solid normal to the growth 
axis was calculated using the formula 



1 



s ( z )= / ^(ip(.x,y,z) + 1) dxdy, 



(10) 



where ip(x,y,z) is the phase-field that varies from +1 in 
the solid to —1 in the liquid in the present model |T^ 1 
and ip(x,y, z) = defines the solid-liquid interface. This 
formula is accurate far enough away from the tip, but of 
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FIG. 2. Plot of the cross-sectional area S(z) vs distance 
\z\ from the tip for different growth parameters. Lengths are 
measured here in units of Wo, the thickness of the diffuse 
interfaces in the phase-field model as defined in Ref. |13|. 



course not at the tip itself, which suffices for the present 
purposes. Plots of S(z) vs \z\ are shown in Fig. || for 
different undercoolings and anisotropy strengths. S(z) 
is accurately fitted by a straight line sufficiently away 
from the tip. This allows us to define the 'parabolic' tip 
radius, p p , of a paraboloid of revolution with the same 
cross-sectional area as the non-axisymmetric phase-field 
shape. Such a paraboloid has a cross-sectional area 



S p (z) — 2irp p (zo — z) for z < Zq, 



(11) 



where zq is its tip position. By fitting the linear part of 
SJz) away from the tip computed from Eq. (|l0|) with Eq. 
(Q l|) , we obtain an accurate estimate of p p (note that zq 
need not coincide with the tip position of the full non- 
axisymmetric shape, see below). The quality of this fit 
will be illustrated below in Fig. 17, where it is compared 



to a different method for determining p p . 

Remarkably, and in good agreement with theoreti- 
cal expectation (Eq. ||), p p coincides to within a nu- 
merical accuracy of 1-2 % with the Ivantsov tip radius 
pi v = 2DPj v /V of the isothermal paraboloid growing 
with the same tip velocity as the phase-field shape for 
the different undercoolings and anisotropics investigated 
here. This indicates that p p w pi v is potentially a good 
scaling parameter for the entire dendrite shape as will be 
confirmed below. Therefore, unless otherwise stated, all 
lengths will be rescaled by pi v for the reminder of this 
paper. The origin of the z-axis is chosen to coincide with 
the tip position of the non-axisymmetric dendrite shape. 

It is also useful to define the equal cross-sectional area 
(ECSA) shape 
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(12) 



where r(z, cf>) is the radial coordinate of the full non- 
axisymmetric shape. This is simply the axisymmetric 
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FIG. 3. Sections of phase-field shapes in the <f> = 0° 
plane (solid line) and <j> = 45° plane (dotted line), equal 
cross-sectional area shape (solid squares), and parabolic fit 
of the latter (dashed line). Parameters are A = 0.1 and 
e 4 = 0.025. 




FIG. 4. Superposition of dendrite contours taken along fins 
and valleys for different parameters. 



shape that has the same cross-sectional area as the full 
shape, and the above results imply at once that r${z) co- 
incides away from the tip with the Ivantsov paraboloid. 
The ECSA shape is shown as solid squares in Fig. p| to - 
gether with the parabolic fit obtained by using Eqs. (|ll|) 
and ( |l^ ) . The parabolic fit is extended all the way to its 
tip position zq for illustrative purposes. 



B. Non-axisymmetric tip morphology 

Longitudinal sections of the dendrite tip in the = 0° 
and <fi = 45° planes that correspond to the 'fins' and 
'valleys' of the non-axisymmetric shape, respectively, are 
superimposed in Fig. ||. The two contours coincide in the 
upper region of the tip, where the shape is essentially ax- 
isymmetric, but depart from a paraboloid (dashed line 
in Fig. ||) in a small region near the tip that will be 
examined more closely in section III.C. The important 
point here, especially relevant for comparisons between 
simulations and experiments, is that aside from this lo- 
calized tip distortion the entire tip shape scales with pi v 
and is ostensibly independent of anisotropy strength and 
undercooling over the range investigated here. This is 
clearly demonstrated by the nearly perfect superposition 
of longitudinal and transverse sections of the various tip 
shapes shown in Fig. || and Fig. || respectively. A 
three-dimensional view of this universal tip morphology 
is shown in Fig. |^. 

In order to test if the analytical form defined by Eq. || 
provides an accurate description of the non-axisymmetric 




FIG. 5. Superposition of transverse tip sections for 
e 4 = 0.025, A = 0.1 (solid lines), e 4 = 0.025, A = 0.2 
(dotted lines), and £4 = 0.0125, A = 0.2 (dashed lines). 
Cross-sections are taken at \z\/pi v = 2, 4, 6, . . .. 
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FIG. 6. Three-dimensional view of the simulated tip mor- 
phology for A = 0.1 and e 4 = 0.025. 

tip shape, we determined A± by minimizing the spatially 
averaged root-mean-square (rms) deviation < 8z 2 > 1//2 , 
between the actual shape in the (OlO)-plane and the 



polynomial form z 



/2 + A^r over the interval 



1 < \z\ < n where n was varied from 4 to 20 in steps of 
2. The resulting values of A\ and < Sz 2 y 1 / 2 are plotted 
for the different undercoolings and anisotropy strengths 
in Fig. i] and Fig. |[ respectively. Fig. || shows that the 
form (|2 ) provides a good fit of the fin shape for n up to 
about 10, after which the rms deviation starts to increase 
rather sharply. Fig. 0, in turn, shows that the fitted value 
of A4 depends on the fitting length. This means that the 
polynomial form (||) does not exactly characterize the fin 
shape. However, the dependence of A 4 on n is rather 
weak: A4 varies between about 0.004 and 0.005 for n be- 
tween 4 and 10. A comparison of the polynomial fit and 
the computed fin shape is shown in Fig. |^, showing an 
excellent overlap up to 10pi v behind the tip. 

We have examined how the form (0) with the 
above procedure to determine A\ fits the entire non- 
axisymmctric shape, and not just the fin. This is il- 
lustrated in Fig. |lC| where we compare the transverse 
sections of the computed phase-field tip shape with the 
ones corresponding to the form (||). The two are in good 
agreement up to ten pj v behind the tip, although a small 
deviation due to higher azimuthal harmonics is notice- 
able. 

We have also fitted the fin shape to the power law 



z = —a\x 



by minimizing the rms deviation from the 
computed shape as above. This minimization was carried 
out both with a and (3 as free parameters and with a as 
free parameter and f3 fixed to the theoretically expected 
value 5/3. The values of a and f3 resulting from these two 
fits are plotted vs the fitting range parameter n in Fig. O 
with the corresponding rms deviations plotted in Fig. ]^ 
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FIG. 7. Plot of A4 vs n obtained from a fourth order poly- 
nomial fit of the fin shape over the interval pi v < \z\ < npi v . 
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FIG. 8. Plot of the root-mean-square shape deviation cor- 
responding to the polynomial and power law fits versus the 
fitting length. Black symbols: £4 = 0.025, A = 0.1; grey 
symbols: £4 = 0.025, A = 0.2; open symbols: £4 = 0.0125, 
A = 0.2. 
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FIG. 9. Comparison of the computed fin shape (open cir- 
cles) with the polynomial fit for n = 10 (A4 = 0.004, solid 
line) and the power law fit for n — 20 and /3 = 5/3 (a = 0.685, 
dashed line). 
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FIG. 10. Comparison of the computed tip cross-sections for 
£4 = 0.025, A = 0.1, and the shape given by Eq. (^|) with 
A4, — 0.004. Cross-sections are taken at \z\/piv = 2, 4, 6, 8, 
and 10. 
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FIG. 11. Plot of the prefactor a and exponent (3 vs the 
fitting range parameter n. Black symbols: £4 = 0.025, 
A = 0.1; grey symbols: £4 = 0.025, A = 0.2; open symbols: 
e 4 = 0.0125, A = 0.2. 



In the two-parameter fit, (3 is larger than 5/3 (ss 1.8) for 
small n and tends to 1.7 for larger n, which is close to 5/3. 
Note that the value of a is somewhat smaller in the two- 
parameter fit than in the one parameter fit because of the 
larger j3 in the former, and that the two-parameter fit has 
a smaller rms deviation as one would expect. The two 
fits, however, become essentially equally good for large 
n. In this range, the one-parameter fit yields a 0.68 
independent of anisotropy strength. A comparison of the 
power law fit and the computed fin shape is shown in 
Fig. 0. In contrast to the polynomial fit, which becomes 
inaccurate for \z\ > 10 pi v , the power law closely fits the 
fin shape even far away from the tip. 



C. Localized tip distortion and tip radius 

Let us now return to examine more closely the shape 
departure from a paraboloid in the region very close to 
the tip. For this purpose, we show in Fig. [l^ a tip 
magnification of the same curves as in Fig. |^, together 
with the axisymmetric shape computed by the boundary 
integral method for an azimuthally averaged surface en- 
ergy. The latter coincides well with the fins and valleys 
in the tip neighborhood, but the superposition of these 
three curves departs from the parabolic fit of the ECSA 
shape, which represents the interface shape if anisotropic 
capillary effects were absent. Note also that the bound- 
ary integral result starts to deviate noticeably from the 
phase-field shape rather close to the tip, approximately 
at \z\ » 0.5 pj v . 

A quantitative measure of this departure is the ratio 
p/piv, where p is the actual numerically computed tip 
radius of the phase-field or boundary integral shape and 
pi v (recall that pi v = 2DP]jJV where Pj v is the Peclet 
number predicted by Eq. [H and V is the numerically 
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computed tip velocity). For the phase- field shape, p is 
computed from the 4> — 0° section using the interpolation 
scheme that is described in Appendix B of Ref. (ll| , and 
which has been tested against exact boundary integral 
benchmark results in two dimensions. 

Figs. [l3| and |lj show the variation p/pi v as a func- 
tion of undercooling and anisotropy, respectively. The 
departure of p from the Ivantsov relation increases with 
anisotropy strength. The same trend was previously 
found for a higher undercooling (A = 0.45) in Ref. |uj 
and the results of these earlier simulations are also shown 
in Fig. |l4|. A new finding here is that this departure be- 
comes independent of undercooling in the low undercool- 
ing range studied here as can be seen from the flattening 
of the curves at small A in Fig. [l3| For A ~ 0.1, p/piv is 
already quasi-independent of A and its small A limiting 
value depends solely on the anisotropy strength. Inter- 
estingly, the variation of p/piv with anisotropy is crudely 
approximated by the relation 



p/piv ~ 1 



(13) 



over an order of magnitude change in anisotropy 
strength, where a — 15e4 is the stiffness anisotropy. 
Since the Gibbs-Thomson condition (^) implies that the 
steady-state dendrite tip temperature is given by 



-2do(l - a)/p, 



(14) 



Eq. |l3| is equivalent to stating that the dendrite tip tem- 
perature at a given undercooling is relatively independent 
of anisotropy strength. The straight line corresponding 
to Eq. |l3| is superimposed as a dashed line in Fig. [l4|. 
The numerical results lie slightly above and below this 
curve for small and large anisotropy, respectively. 
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FIG. 12. Magnification of the tip region showing the same 
curves as in Fig. with the solid squares omitted for clarity. 
Also superimposed is the axisymmetric shape computed by 
the boundary integral method for an azimuthally averaged 
surface energy (thick dashed line). Parameters are A = 0.1 
and e 4 = 0.025. 



IV. DISCUSSION 

A. Comparison with analytical theories 

Let us first compare our results with the analytical 
theories reviewed in the introduction. We have found 
that Eq. U provides a good fit of the tip shape up to a 
distance of about eight to ten pj v behind the tip. More- 
over, for low undercooling, this shape is independent of 
anisotropy strength (at least over the range investigated 
here) which is in qualitative agreement with the predic- 
tion of linear solvability theory (|||. The value of A4 
found here, however, is about twice smaller than the 
value A 4 = 1/96 » 0.104 [§ predicted by this theory. 
One possible reason for this discrepancy is that the ex- 
isting solvability calculations are carried out in the 
limit of vanishing anisotropy, whereas in the present com- 
putations the anisotropy is finite. We have seen, however, 
that even for a short fitting distance behind the tip, A4 
increases from about 0.004 to 0.005 when the anisotropy 
is lowered from 2.5% to 1.25%, which does not appear 



1.0 



0.8 



.0.6 



Q. 



0.4 



0.2 - 



0.0 




O © Boundary integral - e 4 =0.005 

□ E Boundary Integral - £ 4 =0.016 

A A Boundary Integral - £ 4 =0.025 

■ Phase field - e,=0.016 



0.0 



0.2 



0.4 

A 



0.6 
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consistent with an extrapolation of A4 to its theoreti- 
cally predicted value in the limit that €4^0. It seems 
also difficult to explain this discrepancy by the fact that 
existing calculations are based on linearizing the steady- 
state growth equations around a paraboloid of revolu- 
tion. We find indeed a localized tip distortion that de- 
pends strongly on anisotropy strength, and which is not 
accounted for in these theories. The rest of the tip shape, 
however, departs only weakly from a paraboloid and is 
well fitted by the form (Q) independently of anisotropy. 
Therefore, the origin of this discrepancy remains to be 
understood. 

In contrast, the predicted power law (Eq. |j) for the fin 
shape away from the tip is in relatively good quantitative 
agreement with the present simulations. We find a good 
fit of the fin shape with a fixed exponent f3 = 5/3 and 
a prefactor a 0.68 that is independent of anisotropy 
strength and only about 15% larger than theoretically 
predicted, or with a slightly higher exponent (3 = 1.7 
and a lower a « 0.65. It is interesting to note that the 
mapping with two-dimensional growth shapes implies 
that Eq. [| should only strictly hold in a region far behind 
the tip where the cross-sectional shape contains four well- 
developed arms. It is therefore rather remarkable, and 
perhaps coincidental, that it also holds closer to the tip 
in a region where the fins are not yet well developed. 

Finally, we note that our results indicate a smooth 
cross-over from a polynomial form to a power law form 
for the fin shape with increasing \z\ as indicated by the 
excellent overlap of these two forms in an intermediate 
distance range behind the tip (Fig. ph. 



B. Comparison with experiments 

On the experimental side, the most detailed shape 
measurements to date have been carried out by LaCombe 




10 15 20 25 30 35 40 45 



Orientation (j) (deg) 



FIG. 15. Function Q(4>) from the experimental measure- 
ments of LaCombe et al. in SCN |)| (circles and error bars; 
some points close to (j) = 0° have been omitted). The func- 
tion Q{(f>) computed here from the projections of the shape 
defined by Eq. |i| with A4 = 0.004 is shown as a dashed line. 
This curve differs from a pure cosine function, which is also 
shown for comparison. 



et al. who analyzed the full three-dimensional tip mor- 
phology in succinonitrile (SCN) || and by Bisang and 
Bilgram who studied the fin shape in xenon JTo{ ] . 

On the basis of their measurements on SCN dendrites, 
Lacombe et al. reported that the tip shape can be fitted 
by the form 



(15) 



where lengths are in units of the tip radius. The values for 
the function Q(4>) reported by these authors and the as- 
sociated error bars are reproduced here in Fig. [l5|. It can 
be seen that Q(</>) differs from a cosine function, and they 
concluded that the tip morphology cannot be described 
by a pure cos4</> mode, in apparent disagreement with 
solvability theory and the present results. This function, 
however, was constructed from projected shapes obtained 
by taking photographs of the dendrite tip from different 
azimuthal angles </>. We point out here that Eq. 15 with 
Q((j>) obtained from projected shapes does not correctly 
represent the actual non-axisymmetric tip morphology 
from which these projections arc obtained. The reason 
is that the projected shape observed from a given az- 
imuthal angle </> differs from the cross-sectional shape at 
this angle. To illustrate this point, we have drawn in Fig. 
|l6| a cross-section of the tip shape normal to the growth 
direction and rotated it by an angle (f> with respect to 
the viewing direction chosen parallel to the Y"-axis. The 
projected shape appears wider than the 'true' contour of 
the dendrite at that orientation. Therefore, the dendrite 
cross-sections (Fig. 11 in Ref. ||) reconstructed by La- 
Combe et al. using Eq. |l^ are not representative of the 
true cross-sections. 

It is nonetheless possible to relate our results to the 
experimental measurements of these authors by calcu- 
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FIG. 16. Sketch illustrating the difference between the con- 
tour of a dendrite as seen under an angle <p with respect to a 
(OlO)-plane and the 'true' dendrite contour at this angle. 



lating numerically the function Q((p) that corresponds to 
the tip morphology characterized by a single cos 40 mode 
(Eq. §). As mentioned earlier, Eq. | with Ai = 0.004 
provides a good fit to our simulated tip morphology up 
to 10pi v behind the tip and is therefore adequate for the 
purpose of this comparison (neglecting small corrections 
due to higher azimuthal harmonics). To obtain Q (</>), we 
must first calculate numerically the projected shape, and 
thus relate the apparent width w((p, z) of the dendrite 
at a given distance z behind the tip to the true width 
r(— 0, z), which is the intersection of the fixed X-axis 
with the interface. Using the definitions of Fig. and 
purely geometrical considerations, it is simple to obtain 
the relation 

w(4>, z) = r(ip — 4>, z) cos ip (16) 
where ip((p,z) is defined implicitly by the relation 



ip = arctan 



d^rjip - 4>, z) 
r(i/j - (p, z) 



(17) 



Eq. follows directly from the condition that the point 
on the true contour that corresponds to the apparent 
(perceived) edge of the solid tip (point A in Fig. [l6]) be 
the one of maximum width as a function of ip. For a fixed 
0, the corresponding ip can be simply obtained by setting 
equal to zero the derivative of the right-hand-side of Eq. 
|lq with respect to ip. The angle ip defined implicitly by 
Eq. [I?] is then the polar coordinate of this point measured 
from the fixed X-axis in Fig. [If], and (p — ip measures the 
angle between the line OA and the (010) plane. Note that 
at this point, Eqs. 16 and U7\ can be used to construct 



the projected shape that corresponds to an arbitrary true 
shape whose transverse sections are convex. To proceed 
further, we now restrict our attention to the case where 
the true shape is given by Eq. ||, or equivalently 



r((p,z) = 



1 — y/1 — 1QzA± cos' 
4A 4 cos 40 



1/2 



(18) 



Combining Eqs. ^ through we obtain numerically 
the projected shape z(0, w) that is the inverse of w(<p, z) 
at fixed (p. We finally obtain Q((p) by following the same 
procedure as LaCombe et al. that consists of fitting these 
shapes to a fourth order polynomial of the form 



-w 2 /2 - Q(<P)w 4 



(19) 



for different values of (p. This procedure yields a function 
Q((p) that differs from cos 40 and is superimposed as a 
dashed line in Fig. [l5| This function fits well through 
the data points for Q{<p) reported by LaCombe et al. 
We therefore conclude that, within experimental error 
bars, our simulated universal tip morphology dominated 
by a single cos 40 mode with an amplitude A4 w 0.004 
agrees quantitatively well with the true underlying tip 
morphology of SCN dendrites. 

A few additional remarks should be made. Firstly, in 
drawing the above conclusion we have implicitly assumed 
that the inverse problem that consists of reconstructing 
the true shape from its projections has a unique solution. 
It can indeed be shown that this solution is unique p7| , 
such that this assumption is valid. That is if we were to 
calculate the true shape from the projections defined by 
Eq. [l9| we would recover a tip shape that is well fitted 
by Eq. |[ Moreover, an explicit solution to this inverse 
problem can be obtained analytically by cleverly noting 
fll7| that it is exactly analogous to the one of construct- 
ing a two-dimensional equilibrium crystal shape, with the 
role of the anisotropic surface energy being played here by 
w((p, z) at fixed z. This is the well-known Wulff construc- 
tion Jig ] and the application of the £- vector formalism of 
Cahn and Hoffman |19| yields that the distance AB in 
Fig. |l^ is equal to \d^w((p, z)\, and thus that the solution 
to this inverse problem is given by 



r(ip -4>,z) = ^J[w(cP,z)} 2 + [d^w^z)} 2 , (20) 
d4,w(<p, z) 



iP 



-arctan 



w((p,z) 



(21) 



It should therefore be in principle possible to reconstruct 
the higher azimuthal harmonics from the experimental 
data. At present, however, such an analysis is precluded 
by the magnitude of the error bars in the measurements, 
especially in view of the small amplitude of such harmon- 
ics predicted by our simulations. Lastly, we note that in 
the present example the exact (numerically computed) 
projected shapes z((p,w) defined by Eqs. |ll|{l§| are only 
approximately fitted by fourth order polynomials in w, 
even though the underlying true shapes defined by Eq. 
H are exact fourth order polynomials in r. It is simple 
to see why this is so by using Eqs. p|{lS| to derive an 
analytical expression for the projected shape in powers 
of A4 (valid close to the tip), and by noting that this 



expansion generates terms 



at O(Af). Nonetheless, 
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the form (Eq. |T^) provides a reasonably accurate global 
fit of the projected shape over the range < \z\ < 10 and 
is therefore quantitatively adequate to interpret experi- 
mental results. 

In xenon, Bisang and Bilgram have reported that they 
cannot accurately fit the fin shape with a low order poly- 
nomial. They find, instead, a good fit with a power law 
\z\ = ax 13 with a = 0.58 ±0.04 and (3 = 1.67 ±0.05, that 
extends rather close to the tip. Our present results dif- 
fer from theirs in that both for the polynomial and the 
power law fit, the fitting parameters vary with the fitting 
length, including the exponent (3 which increases when 
the fit is restricted to a region close to the tip. On the 
other hand, in agreement with their findings, the poly- 
nomial fit becomes inaccurate for larger fitting lengths, 
whereas the power law fits the fin shape even far be- 
hind the tip. While the exponents obtained from our 
simulations and their measurements are very close, we 
find here a slightly larger prefactor a ~ 0.68, indepen- 
dently of anisotropy. Note, however, that this discrep- 
ancy might be simply due to the choice of the rescaling 
length. Bilgram and Bisang use the actual tip radius as 
a scaling length, whereas we use the Ivantsov tip radius. 
Taking this into account, their prefactor a' and our pref- 
actor a should be related by a' = a{pj pi v Y~ x ■ As p is 
smaller than pi v , our prefactor should indeed be larger 
than theirs. It would be interesting to re-examine the 
experimental data using pi v as the scaling length. 



C. Tip distortion 

The tip distortion analyzed in section III.C can be in- 
terpreted to result from capillary effects at the tip that 
persist even in the limit of low undercooling. It may 
appear at first counter-intuitive that such effects remain 
important in this limit since the magnitude of capillary 
corrections to the tip temperature vanishes due to the 
increase of the tip radius, i.e. d /p — > as A ^ 0. 
One must recall, however, that the tip velocity and thus 
the temperature gradient at the tip also vanishes in this 
limit. Thus the correct measure of the relative impor- 
tance of capillary effects is the ratio of the magnitudes of 
the normal gradient of u in the tip region induced by cap- 
illary variations, (d n u) c ~ do/ ' p 2 , and by heat diffusion, 
(d n u) d ~ V/D, or (d n u) c /(d n u) d ~ Dd /(p 2 V) ~ a*. 
The constancy of a* implies that capillary effects can 
produce a non-negligible tip distortion that persists in 
the low undercooling limit. 

With regard to the tip radius, Eq. [l^ can be crudely in- 
terpreted to be the simplest linear interpolation between 
the isotropic limit where the tip radius approaches its 
isothermal value, and the limit a — * 1 (e^ — > 1/15) which 
marks the appearance of cusps at the [100] orientations 
of the equilibrium shape where p vanishes. There is, of 
course, no obvious reason why this linear interpolation 
should exactly hold in between these two limits and Fig. 



pT| shows that our numerical results do not lie exactly on 
it. Eq. [l3| should only be considered a reasonable first 
estimate of the actual tip radius. 



D. Operating state 

Traditionally, the operating state of the dendrite tip 
has been characterized in terms of two independently 
measurable parameters, the tip radius p and velocity V, 
from which one defines the scaling parameters 



P = 



pV 
2D 
2Dd 
p 2 V 



(22) 
(23) 



It is clear, however, from the present results that an ac- 
curate measurement of the 'true' tip radius is most likely 
not experimentally feasible. Such a measurement would 
require a very high resolution of the localized shape dis- 
tortion near the tip that is already barely noticeable on 
the scale of Fig. [| for a 2.5% anisotropy, and with the 
interface and the reference parabola represented by thin 
lines that are finer than the experimental resolution. 

The present results show that, outside the small region 
very close to the tip where this distortion is noticeable, 
the entire tip shape is well fitted by a paraboloid with a 
small non-axisymmetric four fold deviation, and that the 
rest of the dendrite shape further away from the tip scales 
with the tip radius p p of this paraboloid. Therefore, a 
better definition of the tip operating state is to use p p 
instead of p and to define accordingly the dimensionless 
parameters 



PpV 
2D 
2Dd 

w 



(24) 
(25) 



which are actually the ones that have been traditionally 
measured in experiments. In addition, the definition of 
a* used in solvability theories (see section I) coincides 
with the latter definition with the further assumption 
that p p = pi v . For a weakly anisotropic material such as 
SCN p is not too different from p v such that the two defi- 
nition sets (Eqs. [22]-p3| and p4 -p5|) are roughly equivalent. 
In contrast, for a more strong y anisotropic material such 
as pivalic acid (PVA), p p can be about two to five times 
larger than p, and concomitantly a* four to twenty-five 
times smaller than a* if we assume that £4 is somewhere 
in the range 0.025 — 0.05, where the lower limit has been 
measured by Muschol et al. |20| and the upper one by 
Glicksman and Sin, 
If we adopt Eqs. 



as a definition of the tip operat- 
ing state, three questions remain to be addressed. Firstly, 
what is the best way to measure p p ? Secondly, how ac- 
curately does linear solvability predict the tip operating 
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FIG. 17. Plot of the parabolic tip radius p p versus the fit- 
ting range n (\z\ < npi v ) in a fit using Eq. ^ with both 
p p and Ai as fit parameters. For comparison, the results 
obtained from a fit of the equal cross-sectional area shape 
according to Eq. (fill) are shown as open symbols. 



state as compared to the present simulations ? Lastly, 
does p p necessarily equal pi v , as found here and assumed 
in solvability theory ? 

In our simulations, we have obtained p p in section III. A 
from a plot of the cross-sectional area of the dendrite ver- 
sus distance from the tip, exploiting the fact that, away 
from the tip, p p is simply the slope of this curve divided 
by 2tt. This provides a very accurate procedure, as can 
be demonstrated by plotting the resulting value for p p 
versus the fitting range (Fig. 17): for n between 4 and 
10, p p varies by less than one percent. The dendrite 
cross-section, however, is not easy to measure as the true 
dendrite shape is difficult to reconstruct from longitu- 
dinal projections for the reasons emphasized in section 
IV. B. The best data are obtained for the fin shape |p|,^0|. 
One way to extract p p is therefore to fit the un-rescaled 
fin shape with a fourth order polynomial, 



ir/ Pp ) 2 /2 + A 4 (r/p p ) 4 



(26) 



over a varying distance from the tip where both p p and 
A4 are allowed to vary, which is the method used by 
LaCombe et al. in SCN |J. We have carried out this 
same procedure on our computed fin shapes and the re- 
sults are shown in Fig. [l^ for two different anisotropies. 
(Note that this fit differs from the one carried out in sec- 
tion III.B where we fixed p p to its value (« pi v ) extracted 
from the cross-section measurement, and only varied A 4 .) 
One can see that p p increases with the fitting distance 
(= npiv) behind the tip, which renders a precise determi- 
nation of p p difficult. This trend was observed by Bisang 
and Bilgram in xenon dendrites |10|; on the other hand, 
for SCN, where £4 is about twice smaller than the low- 
est anisotropy studied here, LaCombe et al. find that p p 
is almost constant when n varies between 4 and 10. It 
is therefore possible that this fitting procedure improves 
for lower anisotropies. Fig. 17 also shows that the de- 



TABLE I. Comparison of the selection constants a* ob- 
tained from the present phase-field simulations, from linear 
solvability theory, and from boundary integral calculations. 
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parture of the fitted radius from pi v is of the order of a 



few percent for fitting ranges around n = 10, such that 
this method can provide a reasonable estimate. 

We compare in Table I the values of a* predicted by the 
linear solvability theory of Barbieri and Langer |^3| with 
the values corresponding to the non-axisymmetric shapes 
of the phase-field model and the axisymmetric shapes 
computed by the boundary integral method for an az- 
imuthally averaged surface energy. We list as well the a* 
values corresponding to these shapes. Note that p p = pi v 
both for the paraboloid assumed in solvability theory and 
for the computed non-axisymmetric and axisymmetric 
shapes, such that comparing a* values is equivalent here 
to comparing scaled velocity values Vdo/D = 2a p Pf v . 
One remarkable fact is that the linear solvability the- 
ory predicts relatively accurately a* even though it does 
not describe the localized tip distortion that causes p to 
depart from p p . Thus, we can conclude from this com- 
parison that this distortion does not strongly affect the 
selection of the velocity. A reexamination of the phase- 
field results of Ref. Q in terms of a* show, however, 
that the linear solvability theory becomes increasingly 
inaccurate for larger anisotropy values (£4 > 0.03). 

Finally, McFadden et al. have recently carried out 
a perturbative analysis of the diffusion field around a non- 
axisymmetric isothermal shape defined at leading order 
in A4 by Eq. [2| where A4 is treated as an expansion 
parameter. Since their analysis neglects capillary effects, 
their tip radius p should be compared to p p here. They 
derived a correction to the Ivantsov relation (Eq. |l|) that 
can become significant in the limit of very low undercool- 
ings. For the lowest undercooling studied here (A = 0.1), 
this correction is of the order of a percent and thus com- 
parable to the accuracy at which p p was determined nu- 
merically. Our results are therefore not in contradiction 
with their predictions. We note, however, that the uni- 
versal non-isothermal shape found here starts to deviate 
from the isothermal one they consider only a short dis- 
tance away from the tip where the fins develop under 
the action of anisotropic surface tension. Simulations at 
substantially lower undercooling would be necessary to 
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test if this difference between the isothermal and non- 
isothermal shapes affects their predictions. 



V. CONCLUSIONS 



puting Center and the Northeastern University Advanced 
Scientific Computation Center. We thank Geoffrey Mc- 
Fadden for pointing out to us the analogy between the 
Wulff construction and the reconstruction of the dendrite 
shape from its projections. 



We have studied the three-dimensional morphology of 
the dendrite tip using recent improvements of the phase- 
field method |ll]-|l3| that make it possible to carry out 
quantitatively accurate simulations are relatively low un- 
dercoolings. Our main finding is that the experimen- 
tally measurable low undercooling tip morphology is in- 
dependent of anisotropy strength and thus universal un- 
der the assumption that kinetic effects are negligible, in 
qualitative agreement with solvability theory |3|,|6|. The 
non-axisymmctric deviation of this morphology from a 
paraboloid is well fitted by a single cos 4cj> mode with an 
amplitude that is about twice smaller than predicted by 
solvability theory and in good agreement with existing 
shape measurements in SCN ||. Moreover, these mea- 
surements were reanalyzed here and found to be consis- 
tent with a non-axisymmetric tip shape dominated by a 
single cos4</> mode as in our simulations. 

The fins are well described away from the tip by the 
power law derived by Brener on the basis of the anal- 
ogy @ between three-dimensional steady-state shapes 
and two-dimensional time-dependent growth shapes 
albeit with a slightly larger prefactor than predicted. In- 
terestingly, the validity of this power law extends remark- 
ably close to the tip. Our findings are also in good agree- 
ment with experimental data on xenon dendrites flCfl . 

Finally, we conclude that the 'true' tip radius is not 
an experimentally adequate parameter to characterize 
the tip operating state since the anisotropy-dependent 
shape distortion near the tip that fixes this radius is most 
likely not measurable. In contrast, the tip radius of the 
paraboloid which underlies the rest of the tip morphology 
(excluding this distortion) is both measurable and a good 
scaling parameter for the entire dendrite shape. The lat- 
ter tip radius is indistinguishable from the Ivantsov pre- 
diction over the range of undercooling studied here, which 
does not exclude differences between these two radii to 
be present at even lower undercoolings |p2f . 

It would be interesting in the future to extend the 
present study to investigate how the anisotropic kinet- 
ics of molecular attachment at the interface alters the 
tip morphology. The detailed study of the tip morphol- 
ogy, and in particular its departure from the universal 
shape characterized here, may actually provide a sensi- 
tive probe of interface kinetic effects. 
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